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Introduction 

Composite materials provide design flexibility in that fiber placement and 
orientation can be specified and a variety of material forms and manufacturing processes 
are available. It is possible, therefore, to "tailor" the structure to a high degree in order to 
meet specific design requirements in an optimum manner. Common industrial practices, 
however, have limited the choices designers make. One of the reasons for this is that 
there is a dearth of conceptual/preliminary design analysis tools specifically devoted to 
identifying structural concepts for composite airframe structures. Large scale finite 
element simulations are not suitable for such purposes. 

The present project has been devoted to creating modeling and design analysis 
methodology for use in the tailoring process of aircraft structures. Emphasis has been 
given to creating bend-twist elastic coupling in high aspect ratio wings or other lifting 
surfaces. The direction of our work was in concert with the overall NASA effort Twenty- 
First Century Aircraft Technology (TCAT). A multi-disciplinary team was assembled by 
Dr. Damodar Ambur to work on wing technology, which included our project. 

Summary of Accomplishments 

Our current work has included the following items: 
(1) Analysis to design elastically tailored wings with bend-twist coupling one cross 
section at a time. This work appears in Appendix I. 

Tapered wings may be analyzed with the use of item (1) to discrete spanwise 
wing stations and with the loads known. Our grant monitor, Dr. Damodar Ambur, 
redirected us in 2003 to put a special emphasis on making the analysis of highly 



tapered wings a significant priority and to make it efficient. The result is in items 
(2) and (3). 

(2) Analysis of 3D wings with geometric taper, which includes variation of chord, 
airfoil thickness and cover wall thickness has been created and applied without 
bend twist coupling in the wing. 

(3) Use of simplified design criteria which include the following steps: 

(a) The wing cover thicknesses are defined by spanwise bending strain level. 
This is followed by a check for torsional divergence. 

(b) "Rigid" wing loads are used to predict first order elastic deformations. 

(c) Loading redistributions due to first order elastic deformations are 
evaluated. Effect on bending strain is assessed. 

Use the results of (b) and (c) to set composite ply layups in the wing covers for 
"best" results. This will necessitate subdividing the wing planform into zones of constant 
ply layups. The zones influence the manufacturing process, and, therefore, require 
practical judgment. 

Iteration may be required to define the cover preliminary thickness estimate. 
Items (1) and (2) are complete. Item (3) has been conceived and applied to a fictitious 
Reno Air Racer with no elastic coupling in the wings, but not compared with other 
methodology such as indicated in item (1). 

The differences among (1),(2), and (3) is that (1) requires only structures and 
materials technology. Item (2) and (3) require aerodynamic and mass information. 
Also,We have employed aerodynamic strip theory, which lacks 3D resolution, but which 



should provide reliable preliminary design definition for the structural configuration. 
This work appears in Appendix II. 

AH of our work has been devoted to preliminary design analysis. The objective is 
to give the designer guidance in establishing the structural configuration. Once a viable 
configuration is identified, "fine tuning" to save weight and cost can proceed with the aid 
of large scale numerical simulation and perhaps, in addition, optimization. The key 
element is the definition of the ply layups. 

Concluding Remarks 

The effort that we directed to highly tapered wings at Dr. Ambur's suggestion did 
not permit an evaluation of both taper and bend-twist elastic coupling. A request for 
additional supplemental funding and additional time was submitted to NASA Langley in 
order to reach all objectives. Unfortunately, the TCAT program was cancelled and no 
supplementary funds were provided. 

Two Journal of Aircraft papers were written and presented at AIAA conferences. 
Both papers have been accepted for publication and are under revision. A complete 
record of the project accomplishments appears in Appendices I and II. 
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NOMENCLATURE 

A Torsional stiffness matrix 

a l2 Defined in eq. (14) 

a 26 Defined in eq. (15) 

A;; Laminate membrane stiffnesses for wing box covers, or elements of 

matrix A 

C12 Effective Poisson's ratio, eq. (22) 

Ci6 Twist-camber coupling parameter, eq. (23) 

C„ Beam stiffness matrix 

C44 Torsional Stiffness, eq. (6) 

C45 Coupling Stiffness, eq. (7) 

C 55 Bending Stiffness, eq. (8) 

c A Aerodynamic chord 

c s Chord of wing section structural box 

En Young's modulus of composite material in fiber diretion 

E12 Young's modulus of composite material in transverse direction 

G12 In-plane shear modulus 

H Height of wing section structural box 

h Skin thickness of load bearing covers of the wing structural box 

hk Thickness of k-th ply of a laminate 

Kjj Structural box cover stiffnesses 

Kii Extensional stiffness, eq. (2) 

K12 Shear stiffness, eq. (3) 

iii 



K22 Extension-shear coupling stiffness, eq. (4) 

ky Structural box cover stiffnesses per unit skin thickness 

k<; Camber curvature kinematic matrix defined in eq. (21) 

M x Chordwise bending moment, eq. (17) 

M y Spanwise bending moment, eq. (18) 

N xx Axial running load due to bending in the wing box covers 

N xs Membrane shear flow or stress resultant 

N ss Circumferential (hoop) membrane stress resultant 

Qj: Plane stress stiffnesses for each ply of a laminate 

W xx Bending curvature of wing box structural box, eq. (20) 

x,y,z Cartesian coordinates 

1,2 Transformed coordinates 

P Bend-twist coupling parameter, eq. (10) 

Pi Bend-twist coupling parameter for covers only, eq. ( 1 1 ) 

e Strain in fiber direction 

Exx> £ ss Components of extensional membrane strain in cell wall 

£xi|/, E y¥ Strains in the x^-y^ coordinate system 

Yxs Membrane shear strain in cell wall of wing section structural box 

lie microstrain, strainxlO 6 

<J), X Rate of twist of wing box structural box, eq. (19) 

vj/ Laminate rotation angle with respect to the axis of the structure 
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9 Dominant skin fiber orientation angle 

Superscripts 

i,j,k Refers to indices assuming the values 1, 2, 6 

Superscripts 

fw Refers to front spar web 

k Refers to identifying index 

1 Refers to lower wing box cover 

rw Refers to rear spar web 

u Refers to upper wing box cover 
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INTRODUCTION 

Even many years after the invention of fiber-reinforced composite materials, 
many new uses are being found for them. They are used in many different applications 
because of their high strength-to-weight ratio, stiffness-to-weight ratio and increased 
design flexibility. The goal for designing is usually to reduce weight or cost while 
achieving high performance products. Although composite materials cost more than 
metal materials, they can produce cost competitive structures. 

Elastic tailoring of composite materials utilizes their design flexibility. From 
Rehfield's definition 1 , an appropriate selection of structural concept, fiber orientation, ply 
stacking sequences and blending of materials can achieve a specific design goal. In past 
research, elastic tailoring has been used to influence the aerodynamics of the system 2 . 
This is defined as "aeroelastic tailoring." 

From the definition of Shirk, Hertz and Weisshaar , "Aeroelastic tailoring is the 
embodiment of directional stiffness into an aircraft structural design to control 
aeroelastic deformation, static or dynamic, in such a fashion as to affect the aerodynamic 
and structural performance of the aircraft in a beneficial way. " The example of it is 
Grumman's X-29 technology demonstrator. The X-29's forward swept wing has high 
lift-to-drag ratio and increased resistance to torsional divergence. 

The focus of the present work is mainly on bend-twist coupling in high aspect 
ratio wing boxes. This coupling can increase the stability of the aircraft system and resist 
the tip stall of subsonic wings. Bend-twist coupling is produced by the stiffness of the 
angle plies, which are the off-axis plies oriented at angles with respect to the bending axis 



of the structure. The laminated ply stiffness depends upon the material and fiber 
orientation. The coupling created by the stiffness can create the desired effect. 

Two design strategies are used to understand the bend-twist coupling, which are 
angle ply rotation and laminate rotation. These two design strategies have been applied 
to two types of models. One is the ideal tailored box model from Rehfield's work 3 , in 
which only the upper and lower covers are the load-bearing structural elements. The 
second is similar to the first one but with the additional sides or spar webs being 
modeled. In the latter model, all elements of the box are load-bearing. After extensive 
evaluation, the second model is chosen for our study. The first model tends to under 
estimate torsional stiffness. 

An illustration for a large transport wing section is studied to help understand the 
two design strategies for producing bend-twist coupling. 

DESIGN STRATEGIES 

Preliminary Remarks 

Bend-twist coupling is produced by creating stiffness which is not aligned with 
the bending axis of the primary structure. In wing covers made of laminated composites, 
some off-axis plies will serve to produce the desired effect. This is the context which is 
considered here. 
Angle Ply Rotation 

In this thesis, two design strategies are considered: Angle Ply Rotation and 
Laminate Rotation. The angle ply rotation uses axis-oriented plies and an unbalanced set 
of angle plies. If desired, transverse plies also can be included to increase damage 



tolerance. In this method, bend-twist coupling is produced by the unbalanced layup. The 
angle ply orientation can be varied in order to enhance coupling. As shown in Figure 1, 
the axial plies stay in the same position but the angle of the unbalanced set of angle plies 
varies. The angle, 8, rotates with respect to the axial axis of the structure. 
Laminate Rotation 

The second design strategy is called "Laminate Rotation." This method has been 
applied to the swept forward wing of the X-29 and the analysis of bend-twist coupling in 
Ref. 3. This approach chooses an established layup configuration, which includes axial 
plies and angle plies. Transverse plies also can be included if desired. As shown in 
Figure 2, the entire laminate is rotated from the x-y coordinate system to the x v -y v 
coordinate system with respect to the axial axis of the structure. The coordinates natural 
to the solution of the structure problem are the structure coordinates x, y, whereas the 
principal material coordinates are x ¥ , y w . The rotation angle, \y, is the angle measured 
from the x-axis to the x v -axis. It is recommended to choose a balanced layup 
configuration in order to avoid or minimize warping due to processing. After the 
rotation, the plies are no longer aligned with the axis and bend-twist coupling is 
introduced. 

The angle plies in the Angle Ply Rotation method and the laminate in the 
Laminate Rotation method are rotated until they produce the desired bend-twist coupling. 
It is efficient to have all angle plies in angle ply rotation to produce bend-twist coupling. 
However, there is a disadvantage that arises from the common manufacturing-related 
warping produced by processing thermoset composite materials for unbalanced 



configurations. The advantages of angle ply rotation method over the laminate rotation 
method will inspire innovative utilization. For this method, warping effects due to 
unbalanced configurations will be neglected. 

MODEL AND ANALYSIS 

Model 

A simple model is employed to analyze the wing box using the two design 
strategies. This model is based upon a refinement of the one in Ref. 3, which considers 
load bearing covers and spar webs or sides which are not modeled. The model shown in 
Fig. 3 is used for this study. It factors in the effects of the webs being modeled. As a 
result, this model has contributions from every portion of the box. This approach is more 
realistic. Based upon a consensus of Refs. 4-6, the webs are assumed to be made of 
balanced angle plies. The preferred choice is [±45] plies 4 . In addition, the 
recommended choice for the web wall thickness is less than half of the covers. 

The coordinate system of the box model is shown in Fig. 4, where x and s are 
indicated. Also, the figure shows the ply orientation for the covers. The upper and lower 
covers are tailored identically and asymmetrically. They are mirror images of each other, 
which is the "best" way to create bend-twist coupling. The angle 9 is the angle between 
the dominant off-axis and the ply orientation. The laminates for the upper and lower 
covers can be made of axial plies, angle plies, and transverse plies. Since the s 
coordinate surrounds the box model, the signs for the angles are different for the covers. 
As shown in Fig. 5, the angle for the upper cover is positive, but negative for the lower 
cover. 



Analysis of the Model 

The analysis methodology is based upon the Bemoulli-Euler bending assumption 
applied to Rehfield's theory of thin-walled composite beams'. The warping effects and 
transverse shear deformations are neglected here. The walls of the model carry loads in 
plane stress, and the stress resultants for the walls are denoted Nxx, N xs and N ss . The axial 
running load due to bending is N xx , and N xs is the running load due to torsion (shear 
flow). The circumferential stress resultant, N ss is assumed to be negligible. This implies 
no internal pressure. Consequently, 
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Kn, K12 and K22 corresponds to uniaxial extension, shear, and coupling stiffnesses, 
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The membrane stiffnesses, Ajj, are independent of the stacking sequence, which can be 
determined by simply adding the plane stress stiffnesses, Q { , for each ply. For a 
laminate of N plies, 

A s =i]Qi\ (i, J =1,2,6) (5) 

where h k is the thickness of k l ply. 



Since the webs have contributions of stiffnesses for the present model, the 
expressions of the global stiffnesses are different from the previous model in Ref. 7. The 
global stiffnesses are 
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where C 44 , C 45 and C55 are the torsional, coupling, and binding stiffnesses, respectively. 
In eqs. (6)-(8), the superscripts u and 1 denote the upper and lower covers, respectively, 
and fw and rw denote the front and rear spar webs, respectively. The web wall thickness 
is a fraction of the cover wall thickness, h, which ff is the fraction for front spar web and 
f r is the fraction for the rear spar web. 

Since a balanced configuration is chosen for the webs, the coupling stiffness is not 
affected. Furthermore, the bending effects from the webs are neglected 5 ; the only 
stiffness that is affected is C44, the torsional stiffness, which is increased. 

The stiffnesses require the determination of wall thickness by using layup 
coefficients for a unit thickness in the strain limit equation. It is convenient to define the 
membrane stiffness in terms of thickness 

K.^hk, (i,j = l,2) (9) 

and the bend-twist coupling parameters that are convenient to define 3 
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Design Analysis Methodology for Angle Ply Rotation 

The design criterion used is to limit the axial strain in the fiber direction of a ply. 
Using the limit of fiber strain as a criterion is in use at Northrop Grumman . The wall 
thickness corresponds to one of the ply types reaching the strain limit. Before comparing 
the strains in the fiber direction, the strains in the x-s coordinate system have to be 
determined first. The inversion of Eq. (1) shows the expressions for them. 
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Then the strain transformation equations are used to determine fiber direction 
strains for all ply types. The strain for each ply type with a rotation angle 9 with respect 
to the x-s coordinate system is 
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where a !2 and a 26 are defined as follows 
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The cover wall thickness is determined by comparing the strains for all the ply 
types in the configuration and determine when the limit is reached first. The steps in this 
analysis are indicated in Fig. 6. 
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Design Analysis Methodology for Laminate Rotation 

The stress resultants in the x-s coordinate system are transformed to the x^-y,,, 
coordinate system with an angle vj/ in order to get the strain-stress resultant relationship. 



(16) 



where [A] ' is the inverse of the A matrix. 

. Using the same design criterion as Angle Ply Rotation, the wall thickness also 
corresponds to one of the ply types reaching the strain limit. The strain in the fiber 
directions can be found by using strain transformation equations. They are known by 
transforming the strains in the x v -y v coordinate system to 0° for axial plies, 9 for 
balanced angle plies and 90° for transverse plies. The steps in this analysis are shown in 
Fig. 7. 
Deformation Parameters 

After the wall thickness has been found for the appropriate design strategy, the 
behavior of the structure can be determined. For only torsion and spanwise bending, the 
governing equations are 

m, =<:„♦,,+(:«(- w,„) (17) 
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Eqs. (17) and (18) relate the moments to the rate of twist, <]> iX , and bending curvature, 
W, xx , of the structure. The moments in the x and y directions are M x =2CsHN xs , nose up 
moment, and M y =-N„CsH, spanwise lift (Fig. 3). The rate of twist and bending 
curvature of the structure can be determined from Eqs. (17) and (18) by inversion. 
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where C12 is an effective Poisson's ratio 
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and C26 characterizes the coupling between twist and camber 
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This completes the analysis for each configuration. Note that no special effort has 
been made to tailor for camber response. The focus has been only on bend-twist 
coupling. 



TRANSPORT ILLUSTRATION 

A specific transport illustration is created for analysis using the two design 
strategies. The model shown in Figure 3 with the dimensions of the C-130 aircraft center 
wing box (Fig. 8) is used. The material system for the simple box model is AS4/3501-6 
carbon-epoxy. Its elastic properties are shown in Table 1. The stress resultants in the 
covers are assumed to be 25,000 pounds per inch for N xx , and 5,000 pounds per inch for 
Nxj. The limit for the axial strain in the fiber direction of a ply is set to be 4500 
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microstrain. The configuration for the Angle Ply Rotation design strategy is 50% of axial 
plies and 50% of angle plies; the configuration for the Laminate Rotation design strategy 
is 50% of axial plies and 50% of [±45]. These 50-50 configurations are simplified 
approximations to the wings of the F-18 and the AV-8B aircraft. The transverse plies are 
neglected for both design strategies for simplicity. The webs are fixed at [±45], which is 
based on the unanimous recommendations of Refs. 4-6. For simplicity, the thicknesses 
for the front spar web and the rear spar web are chosen to be the same. The web wall 
thickness is chosen to be forty percent of the cover thickness. The analysis is based on 
this specific wing section and loads. Results are presented for designs using Laminate 
Rotation in Figs. 9-13 and Angle Ply Rotation in Figs. 14-18. 

RESULTS AND DISCUSSION 

Preliminary Remarks 

The illustration corresponds to a transport aircraft wing section, which is that of 
the center wing of the C-130. Thus, the results presented strictly correspond to this 
configuration. This geometry, the chosen material system, design criterion adopted, 
loading considered and strategies used, all contribute to the analytical results and the 
discussion of them. The intention has been to provide an overall, unique assessment of 
wing sections with bend-twist coupling. Additional parametric studies and further 
research are needed to determine if the trends are general. 
Laminate Rotation 

The angle x\f corresponds to the rotation of the laminate with respect to the axis of 
the structure. 
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The following observations are based upon the results in Figs. 9-13. The 
observations are: 

1. The thickness curves have the same geometric shape as the bending and torsional 
stiffness curves with identically placed cusps (Fig. 9 and 10). 

2. The cusps correspond to distinct transitions between ply types reaching the design 
strain limit. If the design limit were raised for the same material system, the 
thickness curve would be lowered, but would maintain the same shape. 

3. If the ratio between N xx and N xs is nearly constant, but with absolute values 
reduced, the thickness curve would be lowered as above. 

4. Changing the mix of plies would likely alter the cusp points, but not the general 
shape of the curves. 

5. Due to cusp placement, the curve of C45 is not exactly skew symmetric (Fig. 10). 

6. The rate of twist and bending curvature curves have the same cusp points but 
opposite curvature for the thickness curve. As a higher thickness corresponds to 
greater stiffness, these trends make physical sense (Fig. 1 1 and 12). 

7. A different trend can be noticed at approximately -60 and -30 degrees in Figs. 1 1 
and 12, respectively. In the neighborhood of -60 degrees, the rate of twist is 
extremely low. In contrast, the bending curvature reaches a minimum near -30 
degrees. 

8. If it is desired to stabilize aeroelastically produced lift due to tailoring, the -60 
degrees is a good choice for Laminate Rotation, since it tends to reduce section 
stall and increase resistance to torsional divergence. 
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9. The camber curvature is small, and it does not follow similar trends to the other 
parameters. Recall that no effort has been made to design for camber (Fig. 13). 

10. More bend-twist coupling does not necessarily correspond to the coupled global 
stiffnesses, as shown in Figs. 11-13. 

Angle Ply Rotation 

The following observations are based upon the results in Figs. 14-20. Due to 
different configurations of ply layups, the results for Angle Ply Rotation cannot be 
directly compared to the results for Laminate Rotation. The angle 9 is the angle ply 
orientation with respect to a fixed [0] axis. The following observations are made: 

1. There is one less cusp to the curves for thickness, stiffness, rate of twist and 
bending curvature. There is a cusp at 0=0°, which could be evaluated by a 
separate analysis. The way the current program is written as shown in Appendix 
A, 0=0° corresponds to dividing by zero, which cannot be determined. 

2. Even though the two design strategies cannot be compared directly, there are still 
some noticeable differences between them. From Fig. 19, the thickness for Angle 
Ply Rotation that is required to satisfy the design cequirements is less than 
thickness for Laminate Rotation. 

3. There is less bend-twist coupling and greater twisting stiffness for negative ply 
rotation (Fig. 15). The bending stiffness is greater for negative angle ply rotation, 
0, than for negative laminate rotation, iy. The curves are affected by the cusps. 

4. Compared to Laminate Rotation, the coupled bending stiffness for Angle Ply 
Rotation is increased and the coupled torsional stiffness is decreased. Based on 
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the changes in stiffnesses, the value of the rate of twist for Angle Ply Rotation is 
greater than Laminate Rotation. 

5. Not only the rate of twist has been affected; the global compliance of bending 
curvature for Angle Ply Rotation is less than Laminate Rotation. 

6. The camber curvature is greater for Angle Ply Rotation than Laminate Rotation 
for positive rotation, but the reverse is true for negative rotations (Fig. 20). 

Concluding Remarks 

An aircraft illustration of the two design strategies has been presented. The trends 
are likely to be affected by selected loading magnitudes and their respective ratios, which 
may be altered by control surfaces and internal features of the structure such as ribs. 

SUMMARY 

In this study, two design strategies specifically chosen to produce spanwise bend- 
twist coupling have been presented and illustrated. The illustration that has been chosen 
is the center wing section of a transport wing. Even though the results of the two design 
strategies show some similarities, the wing cover ply layups are so different that they 
cannot be directly compared. It does appear, however, that the new angle ply rotation 
strategy may produce a lighter weight structure. 

This study is only a beginning. In order to understand these two design strategies 
well, much future work is needed. 
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TABLE 1. NOMINAL UNIDIRECTIONAL MECHANICAL PROPERTIES FOR AS4/3501-6 GRAPHITE-EPOXY 

(ROOM TEMPERATURE, DRY) 

FIBER DIRECTION, TENSION MODULUS, E n 20.0xl0 6 psi 

TRANSVERSE DIRECTION, TENSION MODULUS, E 22 1.7xl0 6 psi 

IN-PLANE SHEAR MODULUS, G 12 0.85xl0 6 psi 

POISSONS'S RATIO, v 12 0.30 
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FIGURE 1 - ANGLE PLY ROTATION 
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FIGURE 2 - LAMINATE ROTATION 
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FIGURE 3 - BOX MODEL WITH WEBS MODELED 
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FIGURE 4 - SIGN CONVENTION FOR PLY ORIENTATION AND LAYUP 
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FIGURE 6 - FLOW DIAGRAM FOR DESIGN ANALYSIS METHODOLOGY 

FOR ANGLE PLY ROTATION 



Rotate the angle plies from the x-s coordinate 
system to the 1-2 coordinate system with an angle, 
0, where #is the angle from the x-axis to the 1-axis. 



Express the strain-stress resultant relationships in 



the 1-2 coordinate system. 



Determine the expressions for the strains in 
the [0] and [0\ plies. 



Limit the axial strain in the fiber direction of a 
ply to 4500 p,£. Compare the strains in the [0] 
and [6\ plies. Determine which one will reach 
the limit strain level first. 



The cover thickness is controlled by the ply 
that reaches the limit strain level first. 



Global stiffnesses and camber curvature 
can then be determined. 
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FIGURE 7 - FLOW DIAGRAM FOR DESIGN ANALYSIS METHODOLOGY 

FOR LAMINATE ROTATION 



Rotate the laminate from the x-s coordinate system 
to the Xyryy, coordinate system with an angle, I//, 
where i//is the angle from the x-axis to the x^axis. 



Transform the fixed stress resultants, Nxx 
and Nxs, from the x-s coordinate system to 
the x^y^ coordinate system. 



i 



Express the strain-stress resultant relationships 
in the x^y^ coordinate system. 



Determine the expressions for the strains 
in the [0], [±6] and [90] plies. 



I 



Limit the axial strain in the fiber direction of 
a ply to 4500 [ie. Compare the strains in the 
[0], [±0] and [90] plies. Determine which 
one will reach the limit strain level first. 



i 



The cover thickness is controlled by the ply 
that reaches the limit strain level first. 



i 



Global stiffnesses and camber curvature 
can then be determined. 
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FIGURE 8 - C-130 BOX DIMENSIONS 
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FIGURE 9 THICKNESS VS. LAMINATE ROTATION 
for (50% [0], 50% [±45], 0.4 web's thickness) 
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FIGURE 10 STIFFNESSES VS. LAMINATE ROTATION 
for (50% [0], 50% [±45], 0.4 web's thickness) 
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FIGURE 11 RATE OF TWIST VS. LAMINATE ROTATION 
for (50% [0], 50% [±45], 0.4 web's thickness) 
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FIGURE 12 BENDING CURVATURE VS. LAMINATE ROTATION 
for (50% [0], 50% [±45], 04 web's thickness) 
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FIGURE 13 CAMBER CURVATURE VS. LAMINATE ROTATION 

for (50%[0], 50% [±45], 0.4 web's thickness) 
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FIGURE 14 THICKNESS VS. ANGLE PLY LAYUP ORIENTATION 
for (50% [0], 50% [8], 0.4 web's thickness) 
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FIGURE 15 BENDING, TWISTING and COUPLING STIFFNESSES VS. ANGLE 

PLY LAYUP ORIENTATION 
for (50% [0], 50% [6], 0.4 web's thickness) 
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FIGURE 16 RATE OF TWIST VS. ANCLE PLY LAYUP ORIENTATION 
for (50% [0], 50% [el, 0.4 web's thickness) 
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FIGURE 17 BENDING CURVATURE VS. ANGLE PLY LAYUP ORIENTATION 

for (50% [0], 50% [9], 0.4 web's thickness) 
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FIGURE 18 CAMBER CURVATURE VS. ANGLE PLY LAYUP ORIENTATION 
for (50% [0], 50% [9], 0.4 web's thickness) 
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FIGURE 19 COMPARATIVE THICKNESS 

for (50% [0], 50% [9], 0.4 web's thickness) 
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FIGURE 20 COMPARATIVE CAMBER CURVATURE 
for (50% [0], 50% [8], 0.4 web's thickness) 
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APPENDIX A 
PROGRAM FOR ANGLE PLY ROTATION 

Variable and Parameter Definitions 

Appendix A is the Fortran 90 program that is used to analyze Angle Ply Rotation. 

The following inputs are required for the program. 

Inputs: 

Width of the box model, Cs 

Height of the box model, H 

Elastic properties of the materials used for the box model, E M , E 22 and G i2 

Limit strain for the design criterion, £ 

Axial loading due to bending, N xx , and due to torsion, N xs 

Configuration of the ply layup - fractions of axial plies, angle plies and transverse plies, 

which are f , fe, f<>o- 
Orientation for the ply layup of the webs, w 
Fraction of the wall thickness for the front spar web, ff 
Fraction of the wall thickness for the rear spar web, fr 

After inputting the values for the above parameters, the program generates the 

following outputs. 

Outputs: 

Bend-twist coupling parameter, p 

Thickness of the covers, h 

Rate of twist, # x: 

Bending curvature, W xx 

Camber curvature, kc 

Global stiffnesses, C44, C55, and C45 

Relative weight, which is compared to a balanced brenchmark reference of [0] and [45] 

plies. 

Notes: All of the inputs are in degrees and English units. 

Program Listing 

c============================== 



c PROGRAM for Angle Ply Rotation 



Program coupling 
real*8Ell,E22,G12,vl2,v21,Qll,Ql2,Q22,Q66,C44,C45,C55 
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real*8 Qbarl lu,Qbarl2u,Qbar22u,Qbarl6u,Qbar26u,Qbar66u,eps 
real*8 Qbarl H,Qbarl21,Qbar221,Qbarl61,Qbar261,Qbar661 
real*8 Qbarl lw,Qbarl2w,Qbar22w,Qbarl6w,Qbar26w,Qbar66w 
real*8 Qbarl lb,Qbarl2b,Qbar22b,Qbar66b 

real*8 Qbl 10,Qbl20,Qb220,Qbl60,Qb260,Qb660 

real*8 Qbl lu,Qbl2u,Qb22u,Qbl6u,Qb26u,Qb66u,Q(4) 

real*8 Qbl 19,Qbl29,Qb229,Qbl69,Qb269,Qb669 
real*8 Qbl H,Qbl2I,Qb221,Qbl61,Qb261,Qb66l 
real*8 Qbl lp,Qbl2p,Qb22p,Qbl6p,Qb26p,Qb66p 
real*8 Qbl ln,Qbl2n,Qb22n,Qbl6n,Qb26n,Qb66n 
real*8 Qbl 15,Qbl25,Qb225,Qbl65,Qb265,Qb665,d,e,betamax 

real*8 kl lu,kl2u,k22u,theta,betasq,Cs,He,sumf,fO,ftheta,r90 

real*8 kl Il,kl21,k221,cl2u,c26u,kcu,thick,thickb,wgain,thetamax 
real*8 kllw,kl2w,k22w,kllb,k22b,betal,hbb,htb,Mx,My,dphi 

real*8thetal,pi,a,c,Nxx,Nxs,al2,a26,al2b,a26b,hb,ht,w,ffw,frw,Wxx 
integer i 
c 

c create the data files for the results « 

c 

open(unit= 1 1 ,file='result.dat',status='unknown') 
open(unit=12,file='resultl.dat',status='unknown') 
open(unit=13,file='result2.dat',status='unknown') 
open(unit=14,file='result3.dat',status- unknown') 
open(unit=15,file='result4.dat',status='unknown') 
open(unit=16,file='result5.dat',status='unknown') 

pi=3.14159 
betamax=0.0 

write(*,*) 'Enter the width for the single beam box in inches' 

read(*,*) Cs 

write(*,*) 'Enter its height in inches' 

read(*,*) He 

write(*,*) "Enter the properties of the material' 

write(*,*)'Ell?' 

read(*,*) El 1 

write(*,*) *E22?* 

read(*,*) E22 

write(*,*) 'G12?' 

read(*,*)G12 

wnte(**)'vl2r 

read(*,*) vl2 

write(*,*) 'strain?' 

read(*,*) eps 

write(*,*) "Enter the axial load Nxx' 

read(*,*) Nxx 

write(*,*) "Enter the shear flow Nxs' 

read(*,*) Nxs 
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write(*,*) "Enter the degree for the angle plies of the webs' 

read(*,*) w 

write(*,*) 'Enter the fraction for the front spar web' 

read(*,*) ff 

write(*,*) 'Enter the fraction for the rear spar web" 

read(*,*) fr 

sumf=0 
c 

c state the fractions for [0] or [theta] or [90] plies in the laminates 
c 

Do while (sumf.NE. 1.0) 

write(*,*) The sum of the fraction for all plies has to be 1' 

write(*,*) 'Enter the fraction for zero degree plies' 

read(*,*) fO 

write(*,*) 'Enter the fraction for theta plies' 

read(*,*) ftheta 

write(*,*) 'Enter the fraction for 90 degrees plies' 

read(V) f90 

sumf=f0+ftheta+f90 
End Do 
c 

c determine the reduced stiffnesses at [0] 
c 

v21=vl2*E22/Ell 

Qll=Ell/(l-vl2*v21) 

Q12=vl2*E22/(l-vl2*v21) 

Q22=E22/(l-vl2*v21) 

Q66=G12 

Q(1)=Q11 
Q(2)=Q12 
Q(3)=Q22 
Q(4)=Q66 

c 

c transform the reduced stiffnesses to different angles 
c 

a=0.0 

call Qmat(a,Q,Qbl 10,Qbl20,Qb220,Qbl60,Qb260,Qb660) 

c=90.0 

call Qmat(c,Q,Qbl 19,Qbl29,Qb229,Qbl69,Qb269,Qb669) 

d=45.0 

call Qmat(d,Q,Qbl 15,Qbl25,Qb225,Qbl65,Qb265,Qb665) 

e=d*pi/180 

Qbarl lb=f0*Qbl 10+ftheta*Qbl 15+f90*Qbl 19 

Qbar 1 2b=f0*Qb 1 20+ftheta*Qb 1 25+f90*Qb 1 29 

Qbar22b=f0*Qb220+ftheta*Qb225+f90*Qb229 

Qbar66b=f0*Qb660+ftheta*Qb665+f90*Qb669 
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kl lb=Qbarl Ib-Qbarl2b**2/Qbar22b 

k22b=Qbar66b 

al2b=(cos(e))**2-Qbarl2b/Qbar22b*(sin(e))**2 

a26b=sin(e)*cos(e) 

hbb=(Nxx/kllb)/eps 

htb=(al2b*Nxx/kl lb+a26b*Nxs/k22b)/eps 

if (hbb.LE.htb) then 

thickb=htb 
else 

thickb=hbb 
end if 
c _ 

c determine the properties for the webs 
c 

call Qmat(w,Q,Qbl lp,Qbl2p,Qb22p,Qbl6p,Qb26p,Qb66p) 
call Qmat(-w,Q,Qbl ln,Qbl2n,Qb22n,Qb!6n,Qb26n,Qb66n) 
Qbarl 1 w=(Qbl lp+Qbl ln)/2 
Qbarl2w=(Qbl2p+Qbl2n)/2 
Qbarl6w=(Qbl6p+Qbl6n)/2 

Qbar22w=(Qb22p+Qb22n)/2 

Qbar26w=(Qb26p+Qb26n)/2 

Qbar66w=(Qb66p+Qb66n)/2 

kl lw=Qbarl Iw-Qbarl2w**2/Qbar22w 

k 1 2w=Qbar 16w-Qbarl 2w*Qbar26w/Qbar22w 

k22w=Qbar66w-Qbar26w**2/Qbar22w 

theta=1.0 
betamax=0.0 
Do i=l,90 

thetal=theta*pi/180 
call Qmat(theta,Q,Qbl lu,Qbl2u,Qb22u,Qb!6u,Qb26u,Qb66u) 
call Qmat(-theta,Q,Qbl H,Qbl21,Qb22l,Qbl61,Qb26l,Qb661) 

Qbarl lu=fO*Qbl 10+ftheta*Qbl lu+f90*Qbl 19 
Qbarl 2u=fO*Qb 1 20+ftheta*Qb 1 2u+f90*Qb 1 29 
Qbar 16u=fO*Qb 160+ftheta*Qb 16u+f90*Qb 1 69 

Qbar22u=f0*Qb220+ftheta*Qb22u+f90*Qb229 
Qbar26u=fO*Qb260+ftheta*Qb26u+f90*Qb269 
Qbar66u=f0*Qb660+ftheta*Qb66u+f90*Qb669 

Qbarl ll=fO*Qbl 10+ftheta*Qbl H+f90*Qbl 19 
Qbar 1 2!=fO*Qb 1 20+ftheta*Qb 1 2l+f90*Qb 1 29 
Qbarl61=fO*Qbl60+ftheta*QM6l+f90*Qbl69 

Qbar221=f0*Qb220+ftheta*Qb221+f90*Qb229 
Qbar26I=fO*Qb260+ftheta*Qb261+f90*Qb269 
Qbar66Nf0*Qb660+ftheta*Qb66l+f90*Qb669 
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kl lu=Qbarllu-Qbarl2u**2/Qbar22u 

k 1 2u=Qbar 1 6u-Qbar 1 2u*Qbar26u/Qbar22u 

k22u=Qbar66u-Qbar26u**2/Qbar22u 

klll=QbarllI-Qbarl21**2/Qbar221 
kl21=Qbarl61-Qbarl21*Qbar261/Qbar221 
k221=Qbar661-Qbar261**2/Qbar221 

betal=kl2u**2/(kl lu*k22u) 

c 

c transform the strains from the xy axis to the angle plies 
c 

al2=(cos(thetal))**2-Qbarl2u/Qbar22u*(sin(thetal))**2 

a26=sin(thetal)*cos(thetal)-Qbar26u/Qbar22u*(sin(thetal))**2 
c 

c determine the thickness 
c 

hb=(Nxx/kl 1 u-Nxs* beta l/kl2u)/(eps*(l -beta 1)) 
ht=((a 1 2/k 1 1 u-a26*beta 1 /k 1 2u)*Nxx+(a26/k22u-a 1 2*beta 1 /k 1 2u) 
& *Nxs)/(eps*(l-betal)) 
c 

c determine the ply that controls the cover thickness 
c 

if(hb.GE.ht)then 

thick=hb 
else 

thick=ht 
end if 
c 

c determine the behavior of the stucture 
c 

Mx=2*Cs*He*Nxs 
My=-Nxx*Cs*He 

cl2u=Qbarl2u/Qbar22u-(Qbar26u*kl2u)/(Qbar22u*k22u) 
c26u=Qbar26u/Qbar22u 

kcu=(2/He)*(cl2u*eps+c26u*Nxs/(thick*k22u)) 
C44=thick*(((k22u+k221)*Cs+k22w*(ff+fr)*He) 
& *(Cs*He/(Cs+He))**2) 

C45=thick*(0.5*(kl2u-kl21)*(Cs*He)**2/(Cs+He)) 
C55=thick*(0.25*(kl lu+kl U)*Cs*He**2) 

betasq=C45**2/(C44*C55) 
if (betamax.LE.betasq) then 

betamax=betasq 

thetamax=theta 
end if 
dphi=(C55*Mx-C45*My)/(C44*C55*(l-betasq)) 
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Wxx=(C45*Mx-C44*My)/(C44*C55*( 1 -betasq)) 
wgain=(2*thick+(ffw+frw)*thick)/(2*thickb+(ffw+frw)*thick)-l 
c 

c write the results to the data files 
c 

write(l 1,*) theta,betasq,betal 
write(12,*)thick,hb,ht 
write(13,*) dphi,kcu,Wxx 
write(14,*) theta,hbb,htb 
write(15,*) thickb,wgain 
write(16,*) C44,C45,C55 
theta=theta+1.0 
End Do 

write(*,*) 'Maximum of beta squared at' 
write(*,*) 'theta=',thetamax 
write(*,*) 'max. beta squared=\ betamax 
End 

* 

c SUBROUTINE for transfroming reduced stiffness matrix 

subroutine Qmat(deg,Q,Qbl l,Qbl2,Qb22,Qbl6,Qb26,Qb66) 

implicit none 

real*8 deg,Qbl l,Qbl2,Qb22,Qbl6,Qb26,Qb66,Q(4) 

real*8 Ql l,Q12,Q22,Q66,degl 

intent(in)::deg,Q 

intent(out)::Qb 1 1 ,Qb 1 2,Qb22,Qb 1 6,Qb26,Qb66 

Q11=Q(1) 

Q12=Q(2) 

Q22=Q(3) 

Q66=Q(4) 

degl=deg*3.14159/180 

Qbll=Qll*(cos(degl))**4+2*(Q12+2*Q66)*(sin(degl))**2 
& *(cos(degl))**2+Q22*(sin(degl))**4 
QbI2=(Qll+Q22-4*Q66)*(sin(degl))**2*(cos(degl))**2 
& +Q12*((sin(degl))**4+{cos(degl))**4) 
Qb22=Qll*(sin(degl))**4+2*(Q12+2*Q66)*(sin(degl))**2 
& *(cos(degl))**2+Q22*(cos(degl))**4 

Qbl6=(Qll-Q12-2*Q66)*sin(degl)*(cos(degl))**3+(Q12-Q22+2*Q66) 
& *(sin(degl))**3*cos(degl) 

Qb26=(Qll-Q12-2*Q66)*(sin(degl))**3*cos(degl)+(Q12-Q22+2*Q66) 
& *sin(degl)*(cos(degl))**3 

Qb66=(Qll+Q22-2*Q12-2*Q66)*(sin(degl))**2*(cos(degl))**2 
& +Q66*((sin(degl))**4+(cos(degl ))**4) 

return 

end subroutine Qmat 
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APPENDIX B 
PROGRAM FOR LAMINATE ROTATION 

Variable and Parameter Definitions 

Appendix B is the Fortran 90 program that is used to analyze Laminate Rotation. 

The following inputs are required for the program. 

Inputs: 

Width of the box model, Cs 

Height of the box model, H 

Elastic properties of the materials used for the box model, En, E 2 2 and G12 

Limit strain for the design criterion, e 

Axial loading due to bending, N xx 

Load due to torsion, N xs 4 

Configuration of the ply layup - fractions of axial plies, balanced angle plies and 

transverse plies. Also, the angle for the balanced angle plies, 0. 
Orientation for the ply layup of the webs, w 
Fraction of the wall thickness for the front spar web, ff 
Fraction of the wall thickness for the rear spar web, fr 

After inputting the values for the above parameters, the program generates the 

following outputs. 

Outputs: 

Bend-twist coupling parameter, p 

Thickness of the covers, h 

Rate of twist, # xx 

Bending curvature, W xx 

Camber curvature, kc 

Global stiffnesses, C44, C55, and C45 

Relative weight, which is compared to a balanced brenchmark reference of [0] and [45] 

plies. 

Notes: All of the inputs are in degrees and English units. 

Program Listing 

c PROGRAM for Laminate Rotation 
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Program coupling 

real*8 Cs,He,El 1, E22,G12 > vl2,theta,betal max 
real*8 sumf,fO,ftheta,f90,v21,Ql l,Q12,Q22,Q66,deglmax 
real*8 a,bn,bp,c,e,betal,C44,C55,C45,kl lu,kl2u,k22u,Q(4) 
real*8 Qbl 10,Qbl20,Qb220,Qbl60,Qb260,Qb660 
real*8 Qbl 19,QbI29,Qb229,Qbl69,Qb269,Qb669 
real*8 Qbl lp,Qbl2p,Qb22p,Qbl6p,Qb26p,Qb66p 
real*8 Qbl ln,Qbl2n,Qb22n,Qbl6n,Qb26n,Qb66n 
real*8 Qb 1 1 1 ,Qbl 2 1 ,Qb22 1 ,Qbl6 1 ,Qb261 ,Qb661 
real*8 Qbl 12,Qbl22,Qb222,Qbl62,Qb262,Qb662 

real*8 Qbarl lw,Qbarl2w,Qbarl6w,Qbar22w,Qbar26w,Qbar66w,ffw,frw 
real*8 w,kl Iw,kl2w,k22w,kl H,kl21,k221,betasq,betasqmax,degsqmax 
real*8 Qfu(3,3),Qfl(3,3),Qbar(3,3),Nxx,Nxs,eps,theta j ^ick^h 1 
real*8 Mx,My,dphi,cl2,c26,cur,kc,Nl l,N22,N12,thickb,wgain 
real*8 h2,h3,theta2,theta3,sl l,s22,sl2,s45p,s45n,detA,AI(3,3) 
integer i 
c 

c create the data files for the results 4 

c 

open (unit=ll,file='resuU.dat , ,status= , unknown') 
open (unit= 1 2,file='result 1 .dat',status— unknown') 
open (unit=13,file= , result2.dat , ,status= , unkno\yn r ) 

open (unit=14,file= , result3.dat',status='unknown') 
open (unit=15,file='result4.dat',status— unknown') 
pi=3.14159 

write(*,*) 'Enter the width for the single beam box in inches' 

read(*,*) Cs 

write(*,*) 'Enter its height in inches' 

read(*,*) He 

write(*,*) T2nter the properties' 

write(*,*)*Ell?' 

read(*,*) El 1 

write(*,*) "E22?' 

read(*,*) E22 

write(*,*)'G12?' 

read(*,*)G12 

write(*,*) VI 2?' 

read(*,*) vl2 

write(*,*) 'strain?' 

read(*,*) eps 

write(*,*) "Enter the axial load, Nxx' 

read(*,*) Nxx 

write(*,*) 'Enter the shear flow, Nxs' 

read(*,*) Nxs 

sumf=0 
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betalmax=0.0 
betasqmax=0.0 
c 

c state the fractions for [0] or [theta] or [90] plies in the laminates 
c 

Do while (sumf.NE. 1.0) 

write(*,*) 'Enter the fraction of zero degree plies' 

read(*,*) fO 

write(*,*) 'Enter the degree for the theta plies' 

read(*,*) theta 

write(*,*) 'Enter the fraction of theta plies' 

read(*,*) ftheta 

write(*,*) 'Enter the fraction of 90 deg plies' 

read(*,*) f90 

write(*,*) "Enter the angle of the webs' 

read(*,*) w 

write(*,*) 'Enter the fraction for the front spar web' 

read(*,*) ff 

write(*,*) 'Enter the fraction for the rear spar web' 

read(*,*) fr 

sumf=f0+ftheta+f90 

End Do 
c 

c determine the reduced stiffnesses at [0] 
c 

v21=vl2*E22/Ell 

Qll=Ell/(l-vl2*v21) 

Q12=vl2*E22/(l-vl2*v21) 

Q22=E22/(l-vl2*v21) 

Q66=G12 

Q(1)=Q11 
Q(2)=Q12 
Q(3)=Q22 
Q(4)=Q66 

c 

c transform the reduced stiffnesses to different angles 
c 

a=0.0 
call Qmat(a,Q,Qbl 10,Qbl20,Qb220,Qbl60,Qb260,Qb660) 

c=90.0 
call Qmat(c,Q,Qbl 19,Qbl29,Qb229,Qbl69,Qb269,Qb669) 

bp=theta 
call Qmat(bp,Q,Qbl lp,Qbl2p,Qb22p,Qbl6p,Qb26p,Qb66p) 

bn=-theta 

call Qmat(bn,Q,Qbl ln,Qbl2n,Qb22n,Qbl6n,Qb26n,Qb66n) 

Qbar(l,l)=f0*Qbl 10+ftheta/2*(Qbl lp+Qbl ln)+f90*Qbl 19 
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Qbar( 1 ,2)=fO*Qb 1 20+ftheta/2*(Qb 1 2p+Qb 1 2n >+f90*Qb 1 29 
Qbar(U)=f0*QbI60+ftheta/2*(Qbl6p+Qbl6n)+f90*Qbl69 

Qbar(2,l)=Qbar(l,2) 
Qbar(2,2)=fO*Qb220+ftheta/2*(Qb22p+Qb22n)+f90*Qb229 
Qbar(2,3)=f0*Qb260+ftheta/2*(Qb26p+Qb26n)+f90*Qb269 

Qbar(3,l)=Qbar(l,3) 

Qbar(3,2)=Qbar(2,3) 

Qbar(3,3)=fO*Qb660+ftheta/2*(Qb66p+Qb66n)+f90*Qb669 
c 

c deterine the properties for the webs 
c 

call Qmat(w,Q,Qbl 1 1,Qbl21,Qb221,Qbl61,Qb261,Qb661) 

call Qmat(-w,Q,Qbl 12,Qbl22,Qb222,Qbl62,Qb262,Qb662) 

Qbarl 1 w=(Qbl 1 1+Qbl 12)/2 

Qbar 1 2w=(Qb 1 2 1+Qb 1 22)/2 

Qbar22 w=(Qb22 1 +Qb222)/2 

Qbarl6w=(Qbl61+Qbl62)/2 

Qbar26w=(Qb26 1 +Qb262)/2 

Qbar66 w=(Qb66 1 +Qb662)/2 

kl lw=Qbarl Iw-Qbarl2w**2/Qbar22w 

kl2w=Qbarl6w-Qbar!2w*Qbar26w/Qbar22w 

k22w=Qbar66w-Qbar26w**2/Qbar22w 

Do i=0,90 
e=i 

thetal=e*3. 14159/180 
theta2=theta*3. 14 1 59/1 80 
theta3=-theta* 3 . 1 4 1 59/ 1 80 
c 

c transform the stiffnesses into a certain rotation angle 
c 

call Tran(e,Qbar,Qfu) 
call Tran(-e,Qbar,Qfl) 
k 1 1 u=Qfu( 1 , 1 ) Qfu( 1 ,2)**2/Qfu(2,2) 

k 1 2u=Qfu( 1 ,3)Qf u( 1 T 2)*Qfu(2,3)/Qfu(2,2) 

k22u=Qfu(3,3)(Qfu(2,3))**2/Qfu(2,2) 

k 1 1 l=Qfl( 1 , 1 )-Qfl( 1 ,2)**2/Qfl(2,2) 

k 1 21=Qfl( 1 ,3)Qfl( 1 ,2)*Qfl(2,3)/Qfl(2,2) 

k221=Qfl(3,3)-(Qfl(2,3))**2/Qfl(2,2) 

c 

c determine the coupling parameter 
c 

betal=kl2u**2/(kl lu*k22u) 

if (beta 1 max. LT. beta 1) then 
beta 1 max =betal 
deglmax=e 
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End If 
c 

c transform the fixed stress resultants to the 1-2 coordinate system 
c 

N 1 1 =cos(theta 1 )**2*Nxx+2*sin(theta 1 )*cos(theta 1 )*Nxs 
N22=sin(thetal)**2*Nxx-2*sin(thetal)*cos(thetal)*Nxs 
N 1 2=-sin(theta 1 )*cos(theta 1 )*Nxx+(cos(theta 1 )* *2 
& -sin(thetal)**2)*Nxs 

det A=Qbar( 1 , 1 )*Qbar(2,2)*Qbar(3,3)+2*Qbar( 1 ,2)*Qbar( 1 ,3) 
& *Qbar(2,3)-Qbar(l,3)**2*Qbar(2,2)-Qbar(2,3)**2*Qbar(l,l) 

& -Qbar(l,2)**2*Qbar(3,3) 

AI(l,l)=(Qbar(2,2)*Qbar(3,3)Qbar(2,3)**2)/detA 

AI(l,2)=(Qbar(l,3)*Qbar(2,3)Qbar(l,2)*Qbar(3,3))/detA 

AI( 1 ,3)=(Qbar( 1 ,2)*Qbar(2,3)Qbar( 1 ,3)*Qbar(2,2))/det A 

AI(2,1)=AI(1,2) 

AI(2,2)=(Qbar( 1 , 1 )*Qbar(3,3)Qbar( 1 ,3)**2)/det A 

AI(2,3)=(Qbar( 1 ,2)*Qbar( 1 ,3)Qbar( 1 , 1 )*Qbar(2,3))/det A 

AI(3,1)=AI(1,3) 

AI(3,2)=AI(2,3) 

AI(3,3)=(Qbar( 1 , 1 )*Qbar(2,2)-Qbar( 1 ,2)**2)/det A 
c 

c compare the strains in different fiber directions 
c 

si 1=AI(1,1)*N1 1+AI(1,2)*N22+AI(1,3)*N12 

s22=AI(2, 1 )*N 1 1 +AI(2,2)*N22+AI(2,3)*N 1 2 

sl2=AI(3,l)*NH+AI(3,2)*N22+AI(3,3)*N12 

hl=sll/eps 

s45p=cos(theta2)**2*sll+sin(theta2)**2*s22+ 
& sin(theta2)*cos(theta2)*s!2 

h2=s45p/eps 

s45n=cos(theta3)**2*sl l+sin(theta3)**2*s22+ - 
& sin(theta3)*cos(theta3)*sl2 

h3=s45n/eps 
c 

c determine the ply that controls the cover thickness 
c j 

if (hl>=h2 .AND. hl>=h3) then 

thick=hl 
else if (h2>=h3 .AND. h2>=hl) then 

thick=h2 
else 

thick=h3 
end if 
if(e.EQ.O)then 

thickb=thick 
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end if 
c 

c determine the behavior of the structure 
c 

wgain=(2*thick+(ffw+frw)*thick)/(2*thickb+(ffw+frw)*thick)-l 
C44=thick*((Cs*He/(Cs+He))**2)*((k22u+k22l)*Cs 
& +(ff+fr)*k22w*He) 

C45=0.5*thick*(Cs*He)**2/(Cs+He)*(kl2u-kl21) 

C55=0.25*thick*Cs*He**2*(kl lu+kl 11) 

betasq=C45**2/(C44*C55) 

Mx=2*Cs*He*Nxs 

My=-Nxx*Cs*He 

dphi=(C55*Mx-C45*My)/(C44*C55*( 1 betasq)) 

cur=(C45*Mx-C44*My)/(C44*C55*( 1 -betasq)) 

c 1 2=Qfu(l ,2)/Qfu(2,2)-(Qfu(2,3)*k 1 2u)/(Qfu(2,2)*k22u) 

c26=Qfu(2,3)/Qfu(2,2) 

kc=(2/He)*(cl2*eps+c26*Nxs/(lhick*k22u)) 

if (betasqmax.LT.betasq) then ■* 

betasqmax=betasq 

degsqmax=e 
end if 
c 

c write the results to the data files 
c 

write(l 1,*) e,betal, betasq 

write(12,*)hl,h2,h3 

write(13,*) e,thick,wgain 

write(14,*) dphi,cur,kc 

write(15,*) C44,C45,C55 
End Do 

write(*,*) 'max. betal' 
write(*,*) deglmax,betalmax 
write(*,*) 'max. betasq* 
write(*,*) degsqmax.betasqmax 
End 



c SUBROUTINE for transforming reduced stiffness matrix 

C ' :=== ===== ======== 

subroutine Qmat(deg,Q,Qbl l,Qbl2,Qb22,Qbl6,Qb26,Qb66) 
implicit none 

real*8 deg,Qbl l,Qbl2,Qb22,Qbl6,Qb26,Qb66,Q(4) 
real*8 Ql l,Q12,Q22,Q66,degl 
intent(in)::deg,Q 

intent(out)::Qbll,Qbl2,Qb22,Qbl6,Qb26,Qb66 
Q11=Q(1) 
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Q12=Q(2) 

Q22=Q(3) 

Q66=Q(4) 

degl=deg*3. 14159/180 

Qbl 1=Q1 l*(cos(degl))**4+2*(Q12+2*Q66)*(sin(degl))**2 
& *(cos(degl))**2+Q22*(sin(degl))**4 
Qbl2=(Qll+Q22-4*Q66)*(sin(degl))**2*(cos(degl))**2 
& +QI2*((sin(degl))**4+(cos(degl))**4) 
Qb22=Q 1 1 *(sin(deg 1 ))**4+2*(Q 1 2+2*Q66)*(sin(deg 1 ))**2 
& *(cos(degl))**2+Q22*(cos(degl))**4 

Qbl6=(Qll-Q12-2*Q66)*sin(degl)*(cos(degl))**3+(Q12-Q22+2*Q66) 
& *(sin(degl))**3*cos(degl) 

Qb26=(Ql l-Q12-2*Q66)*(sin(deg 1 ))**3*cos(degl )+(Q 1 2-Q22+2*Q66) 
& *sin(degl)*(cos(degl))**3 

Qb66=(Qll+Q22-2*Q12-2*Q66)*(sin(degl))**2*(cos(degl))**2 
& +Q66*((sin(degl))**4+(cos(degl))**4) 

return 

end subroutine Qmat '* 



c= 



c SUBROUTINE for tranforming laminate into rotation angle 

subroutine Tran(d 1 ,Qbar,TR2) 
implicit none 

real*8dl,d,Qbar(3,3),TRl(3,3),TR2(3,3) 

intent(in)::dl,Qbar 

intent(out)::TR2 

d=dl*3. 14159/180 

TRl(l,l)=Qbar(l,l)*(cos(d))**2+Qbar(l,2)*(sin(d))**2 

& -Qbar(l,3)*2*sin(d)*cos(d) 

TRl(l,2)=Qbar(U)*(sin(d))**2+Qbar(l,2)*(cos(d))**2 

& +Qbar(l,3)*2*sin(d)*cos(d) 

TR 1( 1 ,3)=Qbar( 1 , 1 )*sin(d)*cos(d)-Qbar( 1 ,2)*sin(d)*cos(d) 

& +Qbar(l,3)*((cos(d))**2-(sin(d))**2) 

TR 1 (2, 1 )=Qbar(2, 1 )*(cos(d))* *2+Qbar(2,2)*(sin(d))* *2 

& -Qbar(2,3)*2*sin(d)*cos(d) 

TRl(2,2)=Qbar(2,l)*(sin(d))**2+Qbar(2,2)*(cos(d))**2 

& +Qbar(2,3)*2*sin(d)*cos(d) 

TR 1 (2,3)=Qbar(2, 1 )*sin(d)*cos(d)-Qbar(2,2)*sin(d)*cos(d) 

& +Qbar(2,3)*((cos(d))**2-(sin(d))**2) 

TRl(3,l)=Qbar(3,l)*(cos(d))**2+Qbar(3,2)*(sin(d))**2 

& -Qbar(3,3)*2*sin(d)*cos(d) 

TRl(3,2)=Qbar(3,l)*(sin(d))**2+Qbar(3,2)*(cos(d))**2 

& +Qbar(3,3)*2*sin(d)*cos(d) 

TRl(3,3)=Qbar(3,l)*sin(d)*cos(d)-Qbar(3,2)*sin(d)*cos(d) 

& +Qbar(3,3)*((cos(d))**2-(sin(d))**2) 
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TR2(l,l)=TRl(l,l)*(cos(d))**2+TRl(2,l)*(sin(d))**2 

& -TRl(3,l)*2*sin(d)*cos(d) 

TR2( 1 ,2)=TR 1 ( 1 ,2)*(cos(d))**2+TR 1 (2,2)*(sin(d))**2 

& -TRl(3,2)*2*sin(d)*cos(d) 

TR2( 1 ,3)=TR 1( 1 ,3)*(cos(d))**2+TR 1 (2,3)*(sin(d))**2 

& -TRl(3,3)*2*sin(d)*cos(d) 

TR2(2,l)=TRl(l,l)*(sin(d))**2+TRl(2,l)*(cos(d))**2 

& +TRl(3,l)*2*sin(d)*cos(d) 

TR2(2,2)=TRl(l,2)*(sin(d))**2+TRl(2,2)*(cos(d))**2 

& +TRl(3,2)*2*sin(d)*cos(d) 

TR2(2,3)=TRl(l,3)*(sin(d))**2+TRl(2,3)*(cos(d))**2 

& +TRl(3,3)*2*sin(d)*cos(d) 

TR2(3,l)=TRl(l,l)*sin(d)*cos(d)-TRl(2,l)*sin(d)*cos(d) 
& +TRl(3,l)*((cos(d))**2-(sin(d))**2) 

TR2(3,2)=TRl(l,2)*sin(d)*cos(d)-TRl(2,2)*stn(d)*cos(d) 
& +TRl(3,2)*((cos(d))**2-(sin(d))**2) 

TR2(3,3)=TR 1 ( 1 ,3)*sin(d)*cos(d)-TR 1 (2,3)*sin(d)*cos(d) 
& +TRl(3,3)*((cos(d))**2-(sin(d))**2^) 

return 
end subroutine Tran 
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Linear Modeling of Aeroelastic Twist of a High Aspect Ratio Wing with Significant 

Taper 

Brett Sikoia, Lawrence W. Rehfield 

Abstract 

The most significant aeroelastic deformation that occurs in optimized cruise flight 
is the elastic twist that occurs about the spanwise direction of the wing. By shifting the 
distribution of the spanwise aerodynamic loads, this aeroelastic twist can have a negative 
impact on induced drag and structural loading. Hence preliminary design tools that can 
characterize this twist, particularly in cases of significant spanwise taper such as 
presented in this paper, are important for beginning the design of efficient aircraft with 
such wings. 

The preliminary design tool presented here makes several simplifying 
assumptions. Airfoil sections are assumed to behave linearly in both an elastic sense and 
a sectional aerodynamic sense. As this is primarily a structural tool, aerodynamic 
behavior is further simplified to assume that individual airfoil sections are assumed to not 
affect each other in the spanwise direction (Trefftz plane). This is classical "strip theory 
aerodynamics" that has been used for high aspect ratio lifting surfaces for generations. 
Furthermore, with composite construction that is balanced with respect to a spanwise 
axis, the axis serves as an elastic axis (EA) which permits bending and torsion to be 
elastically uncoupled. Emphasis is placed here, therefore, on elastic twist. Simple, but 
rational and familiar, design criteria are utilized to size the wing. 



Nomenclature 

EA = Elastic Axis 

c a (r) =Root Aerodynamic Chord 

X =Taper Ratio 

q =Dynamic Pressure 

y ac =Chordwise Separation Distance of the Center of Pressure and EA 

C| =2-D Coefficient of Lift 

ao =2-D Lift Slope of Airfoil 

c m =Constant Coefficient of Moment 

Ciium =2-D Coeffecient of Lift Under Turning Conditions 

f s =Factor of Safety 

kl 1 =Axial Reduced Stiffness per Unit Thickness 

k22 ^Torsional Reduced Stiffness per Unit Thickness 

Introduction 

Aeroelastic response is difficult to assess in preliminary design. Many of the 
tools used to predict aeroelastic response are difficult to use or require large run times 
(such as ANSYS multi-physics engine runs or a NASTRAN run) and hence are limited 
mostly to final analysis rather than true parametric optimization of an aeronautical 
structure. 

It is the intent of this work to present an alternative means of simultaneous 
aerodynamic/structural analysis of high aspect ratio tapered aerodynamic structures 
loaded normal to their spanwise structural axis that can not only give a realistic first 



scaling of the necessary thicknesses of structural box skinning material, but also of its 
subsequent deformation. This formulation is also applicable to swept wings, as only the 
component of the airflow normal to the wing need be considered as causing aerodynamic 
loading, as is commonly done in preliminary swept wing analysis. Further, and most 
importantly, the taper of the wings need not be constant or continuous in any sense, so 
long as a leading edge normal chord can be specified at any point along the span. 

Methodology 

The geometric parameters of the wing (or other lifting surface), and its structural 
composition, must first be described. Secondly, appropriate candidate structural concepts 
and material systems must be identified. 

The material system later illustrated makes use of composite materials, but 
isotropic materials may easily be used. In the case of anisotropic materials, such as 
carbon-composite, reduced stiffnesses must be calculated for laid up plies in a manner 
consistent with sectional analysis'. Also the design strain limit must be specified, a 
common design limit 2 . Finally, the minimum manufacturable wall thickness of the 
material must be established based upon experience. 

The subsequent example and discussion only describes the method for a linear 
taper, but as long as all the parameters can be given as either explicit functions or points 
that may be interpolated, the methodology can be applied. Currently, many wing 
dimensions are reduced dimensions or ratios taken from other dimensions. For the current 
example, the root chord, c a (r), the taper ratio of tip chord divided by root aerodynamic 



chord, X, and the wing half-span. It is assumed that the internal structure scales 
equivalently with the external aerodynamic structure. Consequently, the dimensions of 
the wing structure box must be given at the root, including the thickness of the vertical 
webs. The webs are often sized as a fraction of the thickness of the upper and lower 
covers based upon experience and transverse shear flow. It is assumed that the wing 
structure box is more or less rectangular and a closed section. While the global structure 
is a tapered prism, the individual finite element sections are modeled as untapered. As the 
bending and twisting deformations are elastically uncoupled, aeroelastic twist can be 
estimated independently . 

The structural material is a 50/50 mix of [0] degree spanwise plies and [±45] 
torsionally stiff plies, i.e. [0 2 , ±45], composed of AS4/3501-6 composite material. This 
laminate mimics the F/A-18 and AV-8B without chordwise plies and gives a fairly good 
compromise between bending and twisting stiffness'. 

The dynamic pressure is the only dimensional aerodynamic loading parameter, q. 
The dimensionless aerodynamic loading parameters are the sectional aerodynamic 
parameters: the coefficient of lift at cruise (1-g) loads, c\, the lift slope of the section, ao, 
the perceived maximum multiple of lg that will be experienced in turns, and the angle- 
of-attack independent coefficient of moment, c m . The inclusion of the c m is a departure 
from the earlier work of Reference 4., in addition to several other extensions made here 
of that work. 

The separation of the quarter chord from the elastic center of the structural box is 
given as y ac . Finally the factor of safety, f s , is used to multiply the maximum g-loading to 
provide for reasonable design margin. 



As one last convergence parameter, the number of elements used in the analysis 
must be specified. A consideration of the convergence of the solution is given 
subsequently. 

The analysis of the aerodynamic structure makes use of a series of aerodynamic 
independent sections that interact elastically along the EA. There is one degree of elastic 
freedom per element, twisting about the spanwise axis. Each element experiences a load 
at its outboard tip (node) and twists elastically according to this load. The outboard load 
and the applied aerodynamic load, based on both the initial ci the induced angle of attack 
times ao, are applied at the inboard station (node). The amount of elastic twist 
experienced by each element under the load applied to it is based on its sectional 
torsional stiffness, itself dictated by the structural box and the stiffnesses of the material 
system. In order to remain consistent with this model, the torsional displacement of the 
inboard node (wing root) of the innermost element is fixed at zero (ci stays at its original 
value), while the load applied to the outboard node of the most outboard element is set to 
zero; hence this element sees no relative elastic displacement. The loadings just 
described can also be seen below in Fig.l. 

A system of linear equations results from the equilibrium of the global structure. 
Note that only the twisting moment about the spanwise axis is considered in the 
displacement of the system; while the system does displace in bending parallel with the ci 
loading, the impact on the aerodynamic load distribution is minor, and, hence, it is only 
considered in the context of the bending design strain limit. Once the nodal 
displacements have been found, all of the resultant loadings can be found. 



In the entire system implemented in the coming example, the loadings are actually 
found three times. In the first run, no displacements are allowed and the bending loads 
are found for the worst case scenario of the maximum g-loading multiplied times the 
factor of safety. These are then used in conjunction with the calculated bending stiffness 
to find the skin thickness necessary to reach the bending strain limit at the root 1 ; the 
minimum thickness is also enforced at this step in order to ensure manufacturability. 

In the second run, the loads are calculated in a manner similar to the thickness 
determination loads, only now the Ci is reflective of cruise conditions. These loads are 
used as a baseline for comparison with the loads that are found in the elastically displaced 
case. The final case, of course, finds the loads produced in the wing that has elastically 
deformed to its loads. 

One final interesting calculation that is carried out is the divergence speed 
estimation. By holding constant the dynamic pressure used to calculate the skin thickness 
but varying the dynamic pressure used to calculate the elastic twist of the wing, the 
displacement of the outer tip of the wing can be found as a function of dynamic pressure 
and speed. Thus, the divergence speed can be found graphically by finding where the 
transverse displacement asymptotically approaches infinity. In the following example, 
the result of this calculation is compared to the typical 70-75% chord calculation 5 . 

Example Structure 

The example considered in this paper is the wing of a Reno Unlimited Class Air 
Racer outlined in a preliminary design study and illustrated in Fig. 2. The Reno airport at 
the time of the Reno air races is generally at about nine thousand feet real altitude, and 



competitive racers in the unlimited class generally race at about 644 to 724 km/h (400 to 
450 M.P.H.), which sets the dynamic pressure. The rest of the pertinent data for the racer 
are listed Table 1. Note that the sign convention for the twist is positive nose up. 
Although the taper of this wing is constant, the methodology is still general enough to 
analyze a wing of arbitrary taper, as mentioned earlier. 



Results and Discussion: Reno Air Racer 

1 . As is visible in Fig.3, the combination of the separation of the coefficient of lift 
from the aerodynamic and the negative coefficient moment results in the negative 
coefficient of moment ultimately pushes the induced angles of attack negative. 
Due to aeroelastic effects, the section lift coefficient has decreased by about 
0.018, a decrease of about 15%. This seems unlikely to cause radical problems, 
but there is certainly a perturbation from any initially "rigid" optimized wing 
configuration. Further notice the two inflection points of the twist distribution of 
the wing. The inboard inflection point is mostly the result of the non-constant 
thickness of the section; otherwise the loading tends to scale proportionally to the 
stiffness with the current scaling. The second inflection point is a result of the 
requirement that the rate of change of the twist will reach zero at the outboard 
section, as this section is free of twisting loads. 

2. The twisting load certainly seems affected by the aeroelastic deformation, 
reducing it by about twelve percent, visible in Fig.4. While the twist does not 



directly affect performance, it is interesting to note its decrease and resulting 
impact on the balance of the aircraft. 

3. The bending moment is less affected by the change in sectional angle of attack, 
with a reduction in overall lift in the wing less than 10%, as seen in Fig.5. The 
qualitative curve exhibits no dramatic shift from its original shape. This same 
trend is seen in a relatively minor 1.6" shift in the spanwise center of pressure. 

4. The minimum manufacturability criterion was triggered slightly before the half- 
span mark, visible in Fig. 6. From the sole perspective of rigid bending loads it 
would seem that very small thicknesses would be all that would be required to 
avoid the strain limit in the outboard sections of the wing; however, not including 
the manufacturability requirement produced unrealistic and undesirable outboard 
displacements. It is seen, therefore, that reducing the minimum manafacturable 
thickness would not necessarily save any appreciable weight. 

5. As expected, the torsional and bending stiffness scale with both the span and the 
stiffness, as is evident from the kink at about half-span in Fig. 7. Notice that both 
never reach zero, as this would produce infinite tipwise displacements and loads, 
neither of which is realistic. 

6. As the loading is initially only 2/15 ,h 's that of the turn conditions that set the 
thickness, it is not surprising that the strain remains far below the 4500 
microstrain limit set for the material, visible in Fig. 8. Had the wing been held 
rigid, the strain would have been constant for the initial inboard section. Clearly 
the strains have been lessened due to the induced angles of attack of the sections. 



Note that once the strain was no longer held constant it drops off rapidly; hence 
the strain limit is not reach in the outboard sections. 
7. In a comparison of the current models estimate of the divergence speed with that 
of the 70-75% section model, both predict divergence at flight speeds over 100 
M.P.H. (160 km/h) faster than that of the anticipated flight envelope, as seen in 
fig. 9. By both estimates then it would seem that the current structural 
configuration is a good initial estimate for preventing divergence. 

Convergence Study 

As is common with lower order elements, such as those with only two nodes, a 
large number of nodes are required to achieve satisfactory convergence; in this case of 
the order of 100. Additionally the outermost element is a fairly gross estimation when 
the elements are fairly large, and in fact one can almost draw a line connecting the kinks 
in the station twists that are created by the final rigid element. It would seem logical to 
conclude that holding the last element to a certain much smaller fraction the size of the 
other elements would increase accuracy for the same computational cost. 

Conclusions 

The preliminary design tool presented in this paper can clearly provide useful 
feedback in the preliminary design of elastic high aspect ratio wings of arbitrary taper. 
By this code alone it can be reasonably construed that many higher order models of the 
aerodynamic performance on "rigid" airframes, such as CFD and wind tunnel analysis, 



might fail to fully predict the behavior of the structure, and may be unwarranted without 
including aeroelastic coupled effects. 

From the specific case in question, several other conclusions may be drawn. 
Depending on the magnitude of the coefficient of moment and the coefficient of lift, the 
elastic response of the wing may induce positive or negative changes in the angle of 
attack. The traditional 70-75% span model for divergence prediction is more 
conservative for this illustration. Structural rigidity towards the wingtips is generally not 
as important at preventing failure of the material as keeping deformations and changes in 
load under control. It may also be noted here that by setting the sectional aerodynamic 
coefficients of the outer part of the wing corresponding to control surfaces and to values 
consistent with a control input, control reversal may also be modeled. 

It is the intention of this preliminary design methodology to allow for analysis of 
more diverse planforms than that analyzed here, and, as designed it can be easily applied 
to a host of other lifting surfaces uncoupled elastic bend-twist response. Finally it runs in 
only a few seconds on almost any desktop computer running Matlab. 



Final Remarks 

What is new? 

1 . Analysis efficiency - suitable for preliminary design or conceptual structural 
design on desktop computers 

2. Simple, but rational, established, and validated modeling approach. 
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3. Illustration - extreme taper of wing, realistic parameter values for aircraft type, 
minimum gauge drives design of outboard wing sections (helps reduce tip twist 
and outboard elastic twist effects) 

4. Valid only for wing covers with balanced ply layups. 
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Appendix: Matlab Code 

Source Matlab code for design methodology 
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Tables 



Quantity 


Value 


Units 


q 


19.6*10 J (2.848) 


(Pa)lb/in A 2 


Ca(r) 


1.122(44.19) 


m(inches) 


Wing box width at root 


0.4* c a (r) 


m(inches) 


Wing box height at root 


0.1*c a (r) 


m(inches) 


y ac at the root 


0.2* c a (r) 


m(inches) 


Wing half span 


2.925(115.6) 


m(inches) 


Cl 


0.12 


- 


ao 


5.7 


- 


Cm 


-3.22* 10" 2 


- 


Clturn 


5*c, 


- 


fs 


1.5 


- 


I 


0.39 


- 


kll 


81.5*10 y (ll.*10 6 ) 


Pa(lbf/in 2 ) 


k22 box covers 


21.34*10 6 (3.095*10 6 ) 


Pa(lbf/in 2 ) 


k22 vertical box members 


35.58* 10 6 (5. 16* 10 e ) 


Pa(lbf/in z ) 


vertical box members 
fraction 


0.4 




microstrain criteria 


0.0045 


m/m(in/in) 


minimum thickness 


1.27*10" J (0.05) 


m(inches) 



Table 1 , Design Data for Example Reno Air Racer 
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Figures 
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Figure 1 , Twisting elastic forces on each element 
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Figure 2, Example Reno Air Racer Wing Planform 



14 



Station twist 
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Figure 3, Station Twist 
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Twisting Moment 
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Figure 4, Twisting Moment 
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Figure 5, Bending Moment 
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Skin thickness 
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Figure 6, Skin (Cover) Thickness 
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Figure 7, Stiffnesses 
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x io~ 4 spanwise microstrain, induced cruise conditions 
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Figure 8, Spanwise Microstrain under Cruise Conditions 
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Wingtip Deflection vs. Right speed 
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Figure 9, Torsional Divergence Prediction 
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Con\ergence Study (station twist, degrees) 
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Figure 10, Convergence Study 
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%Brett Sikola 

%NASA funded wing twist approximation 

%initial debug matlab convergence study 

clear all 

%close all 

format short e 

clc 



%first the geometry and the loading must be defined 

%a linear taper will be assumed; dimensions in inches and lbs 

taper^ratio= .39 ; %dimensionless, between 1 and 
root_aero_chord=44 . 190 ; %inches 

root_cover_height=. l*root_aero_chord ; %inches 

root_cover_length= .4*root_aero_chord ; %inches 

root y ac=.2*root aero chord ; %inches 
wing_half_span=115. 6250 ; %inches 
S=wing_half_span* (root_aero_chord+root_aero_chord*taper_ratio) /2 



%in~2 



q=2.848 ; 
ao=5.7 ; 

cm=-3.22e-2; 
cl=.12; 
later) 

cl_turn= 5*. 12 ; 

coordinated) 
fs = 1.5 
es) 



%lb/in"2 

%dimensionless (lift slope) 

%dimensionless (cm, assumed span-wise constant) 
%dimensionless (assumed to be constant, may curve fit ■/ 



%dimensionless, commonly reached value (not ideal ^ 
%dimensionless, factor of safety (1.5 for manned vehicl ^ 



kll=11.82e6 

k22=3.095e6 

k22web=5.16e6 ; 

webf=0.4 

ickness) 

microstrain=4500 / le6 

h_min=.05 

led plies 



%lbf/in~2 (axial stiffness) 

%lbf/in~2 (shear stiffness) 

%lbf/in /v 2 (shear stiffness of the webs alone) 
%dimensionless (web thickness as fraction of cover th if 

%allowable microstrain tnote pure strain conversion) 
%inches, dictated by manafacturability of uncoup \t 



n=100 ; %input ( 'enter the number of elements, 2 or greater'); 
l=wing_half_span/n ; %individual segment length, constant 

a=zeros (n, 1) ; %create empty matrices that the do loop 

b=a; , c=a;, d=a; %and end conditions fill 

y_ac=a; , ca=a; , cs=a; , H=a; , C44=a,-, scale=a; , mom=a; 

%begin specifying the three diagonals of the matrix 

%the b vector is the diagonal entries, 

%while the a and c vetors are the diagonals to the left and right, respectively 

%the d matrix is the remainder matrix 



%coeffecient scaling loop; these will later be used to fill the a,b,c,d vectors 
for i=l:n 

scaled) = 1- (l-taper_ratio) * ( (i-1) /n + 1/2/n) ; %as all linear dimensions scale prop it 
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ortionally 

y_ac (i) =scale(i) *root_^_ac; %they will be scaled using a value 

ca(i) =scale(i) *root_aero_chord; %calculated at each segment 

cs(i) =scale(i) *root_cover_length; 
H(i) =scale(i) *root_cover_height 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%f ind h using minimum manafacturability and turn bending (rigid) %%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

lift_turn(l)=q*ca(l)*l* ( cl_turn*fs ) ,- 
for i=2:n 

lif t_turn(i) =q*ca(i) *1* ( cl_turn*fs ); %rigid segmentary 1 ^ 

if ting 
end 

shear_turn(n) =lif t_turn(n) ; 
for i=(n-l) :-l:l 

shear_turn (i) =lift_turn(i) +shear_turn (i + 1 ) ; %rigid segment rootside shearing moment 
end 

bend_turn(n) =l/2*l*lif t_turn (n) ; 
for i=(n-l) :-l:l 

bend_turn(i) =bend_turn( i+1) +l/2*l*lif t_turn(i) +l*shear_turn (i+1) ; %rigid bending mom ^ 
ent 
end 

for i=l:n 

h(i) =bend_turn ( i ) / (kll*H(i) *cs (i) *microstrain) ; %needed bending stiffness calculation 

if h ( i ) < h_min 

h(i)=h_min ; lenforce minimum thickness 

end 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%vector of thicknesses complete% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%fill out 'd' matrix, torque forces that are constant regardless of aoa 
for i=l:n 

d(i) = (cl*y_ac{i) +cm*ca(i) ) *ca(i) *q*l; 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%% rigid cruise case %%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

mom_rigid(n) =d(i) ■ 
for i=(n-l) :-l:2 

mom_rigid(i) =d(i) + mom_rigid (i + 1) ; % rigid twisting moment 

end 
mom_rigid ( 1 ) =d ( 1 ) +mom_rigid ( 2 ) ,- 

lift_rigid(l)=q*ca(l)*l*cl ; 
for i=2:n 

lif t_rigid(i) =q*ca (i) *l*cl; %rigid segmentary lifting 
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end 

total_lif t_rigid=sum(lift_rigid) ; 

shear_rigid(n) =lif t_rigid(n) ; 
for i=(n-l) :-l:l 

shear_rigid(i) =lif t_rigid(i) +shear_rigid(i + l) ,- %rigid segment rootside shearing mo it 

ment 
end 

bend_rigid(n)=l/2*l*lift_rigid(n) ; 
for i=(n-l) :-l:l 

bend_rigid(i) =bend_rigid(i+l) +l/2*l*lif t_rigid(i) +l*shear_rigid(i+l) ; %rigid bending ^ 
moment 
end 

cp_rigid=bend_rigid ( 1 ) /shear_rigid(l) ; Reenter of pressure, rigid 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% end rigid cruise case %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% find sectional stiffnesses % 

for i=l:n 

C44(i)=cs(i) /v 2*H(i)~2/ (cs(i)+H(i) )"2*(2*h(i) *k22*cs(i) + 2*h (i) *webf *k22web*H(i ) ) ; 

C55(i)=2*h(i) *kll*cs(i) *H(i)"2/4 ; 
end 

% fill out a,b,c vectors % 

for i=3:n 

a(i)=-C44(i-l)/l; %a(l) and a(2) are zero 

end 

b(l)=l; 

for i=2: (n-1) 

b(i)=( C44 (i-1) +C44 (i) ) /l - ao*ca (i) *q*y_ac (i) *1; %note special end conditions 
end 
b(n)=C44(n-l)/l - ao*ca(n) *q*y_ac (n) *1; 

for i=l: (n-1) 

c (i) =-C44 (i) /l; %c(n)=0 as its original value was retained 
end 

%solve a,b,c,d tri-diagonal matrix (i.e. strictly 3 banded sparse matrix) % 
torque_alphas=tridiag(a,b, c, d) ; 

%% Reduce data from the elastic simultaneous case %% 

mom(n)=d(n) + ao*ca(n) *q*y_ac (n) *l*torque_alphas (n) ; 
for i=(n-l) :-l:2 

mom(i)=d(i) + ao*ca (i) *q*y_ac (i) *l*torque_alphas (i) + mom(i+l); % induced twisting m ^ 
oment 
end 
mom ( 1 ) =d ( 1 ) +mom ( 2 ) ; 
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lift(l)=cl*q*l*ca(l) 
for i=2:n 

lif t(i) =q*ca(i) *1* (cl+ao*torque_alphas (i) ) ; % induced segmentary lifting 

end 
total_lift=sum(lif t) ; 

shear (n) =lift(n) ; 
for i=(n-l) :-l:l 

shear (i) =lift (i) +shear (i + 1) ; % induced segment rootside shearing moment 

end 

bend(n) =1/2*1* lift (n) ; 
for i=(n-l) :-l:l 

bend(i) =bend(i+l) +l/2*l*lif t (i) +l*shear (i+1) ; % induced bending moment 
end 

cp=bend(l) /shear (1) ; % induced center of pressure 

%%%%% failure criteria, microstrain and divergence %%%% 
for i=l:n 

strain(i) =bend(i) / (kll*H(i) *cs(i) *h(i) ) ,- %needed bending stiffness calculation 
end 



x_alphas=[0 1/n: 1/n: (1-1/n) 1] ; 
x=l/2/n : 1/n : (1-1/2/n) ; 

q_td=qtd(x,n,C44,y_ac, 1, S,ao) 

%%%%%%%%%%%%% 
%PL0T 0UTPUT% 
%%%%%%%%%%%%% 

figured) %alphas 

torque_degrees=torque_alphas*180/pi ; 

torgue_degrees (n) 

y=[0 ; torque_degrees (2 :n) ; torque_degrees (n) ] ; 

plot (x_alphas,y) ; 

xlabel ( 'dimensionless span');, ylabel (' station twist (degrees)');, title (' Station twist'); 

figure (2) %twisting moment 

plot (x,mom,x,mom_rigid, ' -. ' ) ; 

legend ( ' flexible moment' , 'rigid moment ' ) ; 

title ( "Twisting Moment') 

xlabel ( ' dimensionless span ' ) ; 

ylabel ( 'twisting moment (Ibf-in) ' ) ; 

figure (3) %shear 

plot (x, shear, x, shear_rigid, ' - . ' ) ; 

legend ( ' flexible shear ',' rigid shear') ; 

title ( 'Shear' ) 

xlabel ( ' dimensionless span ' ) ; 

ylabel (' shear (lbf)'); 

figure (4) %bending moment 
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plot (x, bend, x,bend_rigid, '-.' ) ; 

legend) ' flexible bending moment ',' rigid bending moment') ,- 

title ( ' Bending moment ' ) 

xlabel ( ' dimensionless span ' ) ; 

y label ( 'bending moment (lbf-in) ' ) ; 

figure (5) %h, thickness 

plot(x,h) ; 

title ('Skin thickness') 

xlabel ( ' dimensionless span ' ) ; 

ylabel('skin thickness (inches)'); 

figure (6) %C44, C55 

plot(x,C44,x,C55, '-. ' ) ; 

legend ( 'C44' , 'C55' ) ; 

title ( ' stiffnesses ' ) 

xlabel ( ' dimensionless span ' ) ; 

ylabel( 'stiffness (lbf-in~2) ' ) ; 

figure (7) %microstrain 

plot (x, strain) 

title (' spanwise microstrain, flexible case') 

xlabel ( ' dimensionless span ' ) 

ylabel (' strain (in/in)') 

%%%%%%%%% 
%- end -% 
%%%%%%%%% 
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function [result] =tridiag (a, b, c, d) 

%gauss- Jordan tri-diag solver 

%declare an intermediate matrix to which solutions will begin 

%being written. Declare size of matrix to speed execution 

m=length (d) ; 

bb=zeros (m, 1) ; 

cc=bb; 

dd=bb; 

result=bb; 

%the first row is unaffected by the row ops that eliminate the 

%front row in the first loop, so write it in 

bb(l)=b(l) ; 

cc ( 1 ) =c ( 1 ) ,- 

dd ( 1 ) =d ( 1 ) ; 

%the first loop eliminates the front row by taking the row above it 
%and multiplying it by a value that would make it the same as the leading 
%entry in the next row, and then subtracting with this modified row 
for j=2:m 

bb(j)=b(j)-(a(j)/bb(j-l))*cc(j-l) ,• 

cc(j)=c(j) 

dd(j)=d(j)-(a(j)/bb(j-l) )*dd(j-l) ; 
end 

%this leaves one equation in one unknown in the last row, which is solved 

%explicitly 

result (m) =dd(m) /bb(m) ; 

%the row above it then becomes defined, and this is carried through the 

^remainder of the equations, solving the tri -diagonal system 

for j=m-l:-l:l 

result! j ) = (dd( j ) -cc ( j ) *result( j+1) ) /bb(j) ■ 
end 
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function [result] =qtd(x,n, C44 ,y_ac, 1, S, ao) 

%find station at closest to .75 dimensionless span 

found=0; 

for i=l:n 

test=x(i) ; 

if test > 0.70 

if found == 
answeri=i ; 
found = 1 ; 
end 
end 
end 

Kinverse=l/C44{l) ; 
for i=2:answeri 

Kinverse=l/C44 (i) +Kinverse ,- 
end 
K=l/Kinverse ; 

result=K/S/y_ac ( answer i) /ao ; 



